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Abstract 

We  propose  a  new  isogeometric  shell  formulation  that  blends  Kirchhoff-Love  theory  with 
Reissner-Mindlin  theory.  This  enables  us  to  reduce  the  size  of  equation  systems  by  elimi¬ 
nating  rotational  degrees  of  freedom  while  simultaneously  providing  a  general  and  effective 
treatment  of  kinematic  constraints  engendered  by  shell  intersections,  folds,  boundary  con¬ 
ditions,  the  merging  of  NURBS  patches,  etc.  We  illustrate  the  blended  theory’s  performance 
on  a  series  of  test  problems. 

Key  words:  isogeometric  analysis,  NURBS,  shells,  rotation-free,  nonlinear 


1  Introduction 


Reissner-Mindlin  shell  theory,  also  referred  to  as  “thick  shell  theory,”  which  accom¬ 
modates  transverse  shear  deformations,  has  become  the  predominate  theory  used  as 
a  basis  of  finite  element  implementations.  The  main  attribute  of  Reissner-Mindlin 
theory,  as  far  as  finite  element  technology  goes,  is  that  it  is  (^’-conforming,  that 
is,  standard  6'°-continuous  interpolation  functions  are  appropriate  for  representing 
displacements  and  independent  fiber  rotations.  Various  “locking”  phenomena  are 
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a  hindrance  but  are  either  overcome  or  mitigated  sufficiently  by  technologies  that 
were  initially  developed  in  the  1970s  and  refined  by  numerous  researchers  sub¬ 
sequently.  The  basic  four-node  bilinear  Reissner-Mindlin  element,  appropriately 
treated,  has  proved  particularly  efficient  and  robust,  and  is  the  workhorse  of  com¬ 
mercial  finite  element  programs,  such  as  LS-DYNA  [14],  which  is  extensively  uti¬ 
lized  for  automobile  crash  dynamics  and  sheet  metal  forming. 

A  particular  strength  in  this  context  is  the  generality  of  the  formulation  in  that  it  is 
applicable  to  arbitrarily  large  deformations  and  fully  nonlinear  elastic  and  inelastic 
constitutive  behavior.  It  is  fair  to  say  that  Reissner-Mindlin  elements  have  almost 
completely  dominated  finite  element  shell  analysis  for  over  30  years. 

Historically,  initial  efforts  to  develop  shell  elements  were  based  on  Kirchhoff-Love 
theory,  also  known  as  “thin  shell  theory.”  Kirchhoff-Love  theory  invokes  the  Kirch- 
hoff  constraint  of  zero  transverse  shear  strains,  and  requires  C1  -continuous  finite 
element  basis  functions  for  the  transverse  displacements.  Some  functions  of  this 
type  exist,  but  these  lead  to  extremely  complex  elements  that  are  difficult  to  use  in 
a  general  nonlinear  computing  environment.  Ultimately,  these  elements  could  not 
compete  with  the  far  simpler,  more  efficient  and  robust  Reissner-Mindlin  elements. 
Within  the  traditional  finite  element  paradigm,  C1  -continuous  basis  functions  re¬ 
quire  derivative  degrees  of  freedom. 

With  the  advent  of  isogeometric  analysis  [19,12]  a  plethora  of  C1  -continuous  basis 
functions  is  available  and  new  ones  are  under  development  (see  e.g.  [11,27,25]). 
These  offer  unique  finite  element  possibilities  in  that  their  smoothness  properties 
do  not  involve  derivative  degrees  of  freedom.  Consequently,  so-called  “rotation- 
free”  elements  become  a  possibility. 

This  idea  was  first  proposed  in  Hughes  et  al.  [19]  and  has  been  developed  by  Kiendl 
et  al.  [21,20]  and  Benson  et  al.  [6].  As  the  name  implies,  rotation-free  Kirchhoff- 
Love  elements  do  not  involve  rotational  degrees  of  freedom.  This  represents  a  re¬ 
duction  in  the  total  number  of  degrees  of  freedom  and  also  the  bandwidth  size, 
by  roughly  a  factor  of  two  each.  In  implicit  analysis,  compared  with  the  use  of 
Reissner-Mindlin  elements,  storage  of  the  left-hand-side  jacobian  matrix  will  be 
reduced  by  approximately  a  factor  of  four,  and  factorization  cost  will  be  reduced 
by  approximately  a  factor  of  eight.  These  savings  are  very  significant  and  poten¬ 
tially  could  represent  a  decisive  advantage  when  compared  to  the  use  of  Reissner- 
Mindlin  elements.  In  explicit  analysis,  storage  would  be  reduced  somewhat,  but 
not  by  a  significant  enough  amount  to  be  advantageous,  and  compute  cost  might  be 
commensurate,  but  slightly  greater  than  Reissner-Mindlin  elements,  due  to  the  ne¬ 
cessity  of  computing  second-derivatives,  or  equivalents  (see  Benson  et  al.  [6]  for  a 
variant  not  involving  directly  computing  second  derivatives).  Nevertheless,  the  po¬ 
tential  advantage  in  implicit  analysis  is  so  significant  that  it  warrants  further  inves¬ 
tigation.  Unlike  the  case  for  Reissner-Mindlin  elements,  areas  in  the  shell  that  are 
only  geometrically  C°-continuous  require  special  treatment.  Physically,  these  are 
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associated  with  sharp  folds,  creases  and  stiffener  junctions.  Along  these  C°  lines, 
the  angle  between  intersecting  branches  needs  to  be  maintained  during  deforma¬ 
tion.  Various  ways  of  introducing  these  constraints  have  been  investigated.  Kiendl 
et  al.  [20,21]  have  investigated  different  implementations  of  linear  constraint  equa¬ 
tions  and  penalized  “bending  strips.”  There  are  other  possibilities  that  heretofore 
have  not  been  investigated,  namely,  “rigid  bodies”  [7].  Another  need  for  the  use  of 
these  techniques  is  in  multi-patch  NURBS  meshes  of  even  smooth  shells  in  which 
basis  functions  at  the  patch  boundaries  only  possess  U°-continuity. 

It  is  our  opinion  that  none  of  these  constraint  techniques  is  completely  satisfactory 
in  all  cases.  Even  if  they  can  be  used  successfully  in  certain  situations  they  can  give 
rise  to  reduced  critical  time  steps  in  explicit  dynamic  analysis  and  preclude  plastic 
hinge  formation  at  shell  intersections  and  at  boundaries,  where  often  there  are  very 
large  bending  moments.  The  plastic  hinges  are  precluded  by  the  fact  that  there  are 
no  transverse  shear  strains  in  rotation-free  elements.  As  a  result,  in  situations  such 
as  these  Reissner-Mindlin  elements  are  much  more  effective,  but  these  of  course 
require  rotational  degree  of  freedom.  Our  conclusion  is  that  the  potential  advan¬ 
tages  of  rotation-free  elements  in  implicit  analysis  are  presently  difficult  to  realize 
in  most  engineering  situations. 

In  this  work  we  offer  what  we  believe  is  a  novel  and  practical  solution  to  these 
difficulties.  It  is  based  on  the  concept  of  a  “blended  shell  theory”  in  which  the  as¬ 
sumptions  underlying  classical  Reissner-Mindlin  and  Kirchhoff-Love  theories  are 
combined.  Specifically,  we  describe  a  blended  shell  through  a  linear  combination 
of  the  kinematical  assumptions  of  Reissner-Mindlin  and  Kirchhoff-Love  theories. 
We  recall  that  the  only  difference  between  the  two  theories  involves  the  description 
of  the  kinematics  of  the  fiber  vectors.  In  the  case  of  Kirchhoff-Love  theory,  the  fiber 
vector  remains  orthogonal  to  the  lamina  reference  surface  throughout  deformation 
(see  Hughes  and  Liu  [17]  for  a  description  of  the  concepts  of  shell  laminae  and 
fibers),  whereas  in  Reissner-Mindlin  theory  it  depends  on  the  independent  rotation 
fields. 

The  parameter  defining  the  linear  combination  is  a  function  of  the  coordinates,  but  it 
is  convenient  to  implement  it  discretely  in  terms  of  control  variables  (i.e.,  “nodes”). 
Our  blended  theory  can  represent  both  extremes,  that  is,  an  entire  Reissner-Mindlin 
shell,  or  an  entire  Kirchhoff-Love  shell.  It  can  also  represent  a  judicious  combi¬ 
nation  of  the  two  to  deal  with  the  issue  of  constraints.  In  this  case,  we  propose 
to  select  all  the  control  variables  where  there  are  intersections,  boundaries,  folds, 
or  C°  interfaces  as  Reissner-Mindlin  control  variables,  and  elsewhere  in  the  shell 
as  Kirchhoff-Love  variables.  The  result  provides  for  the  elimination  of  rotational 
degrees  of  freedom  in  the  smooth  regions  of  the  shell  and  the  robust  treatment  of  in¬ 
tersections,  boundaries,  folds,  etc.,  provided  by  Reissner-Mindlin  theory  precisely 
where  it  is  needed.  We  believe  this  represents  an  ideal  solution  that  may  provide 
optimally  efficient  shell  discretizations  for  a  wide  range  of  practical  problems.  An 
application  in  sheet  metal  forming  is  schematically  illustrated  in  Ligure  1 . 
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for  sharp  corners. 


Fig.  1.  Schematic  illustration  of  an  application  of  the  blended  shell  theory  for  an  implicit 
simulation  of  metal  forming.  Regions  with  low  curvature  may  be  efficiently  analyzed  with 
rotation-free  thin  shells  while  regions  of  high  curvature  would  use  the  Reissner-Mindlin 
formulation. 

The  idea  of  a  blended  shell  theory  can  go  even  further.  Specifically,  higher-order 
shell  theories  involve  more  complex  descriptions  of  fiber  behavior.  These  can  be 
blended  with  Kirchhoff-Love  and/or  Reissner-Mindlin  theories  to  achieve  efficient 
formulations  for  composite  laminates  and  sandwich  shells,  and  they  may  also  be 
useful  in  biomechanical  modeling  (e.g.  arteries).  It  is  also  possible  to  go  still  fur¬ 
ther  with  the  concept  and  blend  shell  theories  with  continuum  solid  theories  in  the 
spirit  of  transition  elements.  Again  biomechanical  modeling  opportunities  present 
themselves,  such  as  for  heart-artery  models.  We  also  note  that  all  these  blended 
theories  can  be  developed  within  the  IGA  format  of  exact  CAD  modeling. 

The  blended  formulation  presented  here  is  valid  for  a  broad  class  of  basis  functions 
with  C 1  or  greater  continuity,  but  in  our  examples  it  is  applied  to  B-splines  and 
NURBS.  Their  form  is  summarized  in  Section  2.  This  is  followed  in  Section  3 
by  a  brief  summary  of  current  methods  for  joining  patches  in  rotation-free  shell 
analysis.  A  detailed  presentation  of  how  the  blended  formulation  may  be  efficiently 
implemented  is  provided  in  Section  4.  Linear  and  nonlinear  example  calculations 
follow  in  Section  5  to  show  the  strengths  and  weaknesses  of  the  various  approaches 
to  multi-patch  analysis,  and  the  good  behavior  of  the  blended  formulation  in  all 
cases. 


2  A  summary  of  B-spline  and  NURBS  basis  functions 


The  literature  on  the  basis  functions  used  in  CAD  and  CG  is  extensive.  The  standard 
reference  for  NURBS  is  by  Piegl  and  Tiller  [23],  and  a  popular  introductory  book 
is  by  Rogers  [24],  The  first  paper  on  isogeometric  analysis  [19]  contains  a  concise 
introduction  to  B-splines  and  NURBS,  and  the  definitive  treatment  of  isogeometric 
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analysis  is  provided  by  Cottrell,  Hughes  and  Bazilevs  [12].  A  brief  summary  is 
provided  here  to  establish  the  notation  used  in  the  remainder  of  the  paper. 


A  B-spline  curve  is  defined  by  a  knot  vector  and  a  vector  of  control  points.  The  knot 
vector  defines  the  set  of  basis  functions  N\(s)  in  terms  of  the  parametric  coordinate 
s  through  the  recursive  relation 

Mb  ,  s  I  lif  A4<s<  A4+1 

NaAs)  =  \  (!) 

I  0  otherwise. 


8  SA  xlp-M  +  SAW  - 

SA+P  —  Sa  —  Sa+1 


(2) 


where  p  >  0  is  the  degree  of  the  B-spline.  For  brevity,  the  degree  is  omitted  in  the 
remainder  of  the  paper. 


The  knot  vectors  used  here  are  all  open.  The  first  and  last  knots  have  multiplicity  p+ 
1  for  a  B-spline  of  polynomial  degree  p.  Knots  may  be  repeated  in  the  interior  of  the 
knot  vector,  with  each  repetition  locally  lowering  the  degree  of  continuity  by  one. 
The  locations  of  the  knots  define  the  boundaries  of  the  elements  in  the  parametric 
space.  Note  that  the  basis  functions  are  in  general  not  interpolatory  except  at  the 
boundaries,  but  they  do  satisfy  the  partition  of  unity  property,  J2a  ^a(s)  =  1,  and 
are  non-negative  everywhere. 


A  basis  function  of  degree  p  spans  up  to  p  +  1  elements.  Starting  from  the  left 
boundary,  the  first  basis  function  spans  one  element,  the  second  spans  two  elements 
and  so  on  until  the  basis  functions  span  p  +  1  elements.  The  same  progression 
of  spanning  from  1  to  p  +  1  elements  occurs  working  to  the  left  from  the  right 
boundary. 

The  control  point  xa  is  associated  with  the  basis  function  N\  to  define  the  curve 

x(s)  =  J2na(s)xa-  (3) 

A 

Note  that  the  control  points  do  not  necessarily  lie  on  the  curve  they  define,  except 
for  the  first  and  last  ones. 


Surfaces  are  defined  by  a  Cartesian  product  of  one-dimensional  B-spline  curves, 

x(s1,s2)  =  ^N^b(si,s2)xab  where  s2)  =  A^(si)A^(s2).  (4) 

AB 

For  simplicity,  the  subscript  AB  is  replaced  by  a  single  subscript  in  the  remainder 
of  the  paper.  As  a  result  of  their  construction,  the  two-dimensional  basis  functions 
are  in  general  not  interpolatory. 
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A  NURBS  curve  is  defined  by  introducing  a  fourth,  homogeneous  coordinate  w  > 
0  at  each  control  point. 


N\(s)wa 
■  N^(s)wB/ 


(5) 


The  fourth  coordinate  is  not  an  unknown,  but  specified  with  the  initial  coordinates 
Xa  to  define  the  reference  configuration,  and  it  remains  constant  during  the  analy¬ 
sis.  By  inspection,  the  NURBS  basis  functions  N A  are 


Na(s) 


N^(s)wa 

Ebnb(s)wb' 


(6) 


NURBS  surfaces  are  again  defined  by  a  Cartesian  product  of  the  one-dimensional 
basis  functions,  and,  again,  the  double  subscript  is  replaced  by  a  single  subscript, 

x(s1,s2)  =  J2Na(si,s2)xA-  (7) 

A 

NURBS  curves  and  surfaces  have  the  same  properties  as  B-spline  curves  and  sur¬ 
faces: 

(1)  They  form  a  partition  of  unity,  J2A  NA  =  1. 

(2)  The  support  of  each  NA  is  compact,  spanning  up  to  p  +1  knot  intervals  (ele¬ 
ments). 

(3)  The  basis  functions  are  non-negative.  NA  >  0. 

(4)  The  basis  functions  are  not  usually  interpolatory. 

(5)  The  interior  of  the  patch  may  be  continuous  up  to  Cp~ 1  but  the  continuity 
between  patches  is  C°  unless  constraints  are  imposed  to  increase  it. 

(6)  The  coordinates  of  the  control  points  do  not  necessarily  lie  on  the  surface. 


3  Multi-patch  analysis 


NURBS  and  B-spline  patches  are  naturally  assembled  in  a  manner  similar  to  tradi¬ 
tional  finite  elements  to  build  domains  that  are  more  topologically  complicated  than 
a  logically  rectangular  mesh.  While  NURBS  and  B-splines  can  be  up  to  Cp~ 1  on 
their  interiors,  they  meet  at  C°  boundaries.  For  the  analysis  of  solids,  or  for  shear 
deformable  structural  formulations,  this  does  not  pose  any  problem  because  their 
formulation  requires  only  C°  continuity.  Thin  shell  formulations,  on  the  other  hand, 
require  C1  continuity  to  transmit  bending  moments,  a  condition  which  is  violated  at 
the  patch  boundaries.  The  C°  boundaries  between  patches  behave  like  piano  hinges 
for  thin  shell  formulations,  and  therefore  C 1  continuity  must  be  enforced  along  the 
shared  boundary  to  transmit  the  moments. 
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Joining  patches  edge-to-edge  to  form  a  single  C 1  surface  is  the  simplest  type  of  C° 
merge  occurring  in  structural  analysis.  Much  more  common  are  the  cases  where 
the  final  structure  retains  the  C°  lines,  for  example,  as  in  the  square  tube  buckling 
problem  shown  in  Section  5.4.  In  the  case  of  the  square  tube,  there  are  sharp  90° 
comers  that  must  be  maintained,  making  the  C°  edge  line  a  genuine  geometric 
feature.  A  simple  T-intersection  is  another  case  of  a  C°  line  that  is  also  frequently 
encountered  in  real  structures.  Like  a  sharp  comer,  the  low  level  of  continuity  is 
an  intrinsic  property  of  the  exact  geometry  that  cannot  be  eliminated  by  a  different 
choice  of  basis  functions.  The  problem  for  merging  NURBS  patches  that  define  a 
smooth  surface  can  be  addressed  by  using  different  basis  functions,  e.g.,  T-Splines 
[27],  but  the  problems  with  sharp  edges  and  corners  remain. 

There  are  a  number  of  different  ways  to  address  these  issues.  For  example,  the 
rotation-free  formulation  by  Onate  and  Flores  [22]  uses  6'°  basis  functions  in  a 
novel  way.  Within  the  context  of  isogeometric  analysis,  linear  constraints  formu¬ 
lated  in  a  master-slave  form  [21]  have  been  used.  The  linearity  stems  from  the 
assumption  that  the  membrane  deformations  and  rigid  body  rotations  are  small, 
conditions  that  are  not  satisfied  in  general  for  nonlinear  problems. 

Rigid  bodies  [7]  appear  to  be  the  natural  generalization  of  linear  constraints  to  the 
nonlinear  regime,  offering  a  symmetry  in  their  treatment  of  adjacent  patches  that 
has  been  found  useful  in  contact.  The  coordinates  of  all  points  in  a  rigid  body  are 
constant  in  the  body’s  local  reference  frame.  If  the  points  satisfy  a  linear  constraint 
in  the  local  reference  frame  at  t  =  0,  then  the  constraint  will  be  satisfied  at  all  later 
times  in  the  local  frame  independent  of  any  large  rigid  body  rotations.  Instead  of 
spending  a  lot  of  time  deriving  potentially  complicated  linear  constraints  between 
adjacent  patches  (e.g.,  for  the  L-shaped  domain  used  in  an  Section  5.3),  it  is  suffi¬ 
cient  to  simply  group  each  independent  set  of  control  points  into  a  rigid  body.  For 
example,  if  two  patches  share  ten  control  points  along  a  common  boundary,  there 
are  ten  sets  of  three  linearized  master-slave  constraints,  each  having  three  control 
points  (the  shared  control  point  plus  an  additional  one  from  each  patch).  The  ten 
sets  are  used  to  define  ten  rigid  bodies,  each  consisting  of  the  appropriate  set  of 
three  control  points.  The  particular  rigid  body  formulation  used  in  the  example  cal¬ 
culations  is  summarized  in  the  Appendix. 

Care  must  be  exercised  in  introducing  constraints.  As  shown  in  Section  5.2,  even 
simple  eigenvalue  problems  may  have  significant  sensitivity  to  the  choice  of  con¬ 
straint  equations.  Generating  the  correct  constraints  for  general  structural  topolo¬ 
gies  and  basis  functions  in  an  automated  manner  is  anticipated  to  be  a  difficult 
problem. 

Another  approach  is  the  penalty  method,  e.g.,  as  embodied  in  the  bending  strip 
method  [2].  It  is  assumed  that  the  shell  structure  is  comprised  of  smooth  subdo¬ 
mains,  such  as  NURBS  patches,  that  are  joined  with  G'°-continuity.  In  addition,  thin 
strips  of  fictitious  material,  also  modeled  as  surface  NURBS  patches,  are  placed  at 


7 


structural  patch  intersections.  The  triples  of  control  points  at  the  patch  interface, 
consisting  of  a  shared  control  point  and  one  on  each  side,  are  extracted  and  used  as 
a  control  mesh  for  the  bending  strips.  The  parametric  domain  of  each  bending  strip 
consists  of  one  quadratic  element  in  the  direction  transverse  to  the  interface  and, 
for  simplicity  and  computational  efficiency,  as  many  linear  elements  as  necessary 
to  accommodate  all  the  control  points  along  the  length  of  the  strip.  The  material  is 
assumed  to  have  zero  mass,  zero  membrane  stiffness,  and  non-zero  bending  stiff¬ 
ness  only  in  the  direction  transverse  to  the  interface.  Bending  strips  transmit  the 
bending  moment  from  one  patch  to  the  other,  which  may  not  be  accomplished  with 
a  C°-continuous  discretization  alone.  As  such,  the  bending-strip  method  is  viewed 
as  a  physically-motivated  penalty  technique.  The  bending-strip  method  was  suc¬ 
cessfully  employed  in  the  modeling  of  wind  turbine  blades  in  [2,1]. 

Constraints  are  typically  enforced  by  formulating  them  in  a  master-slave  formalism 
that  imposes  the  constraints  sequentially  when  explicit  time  integration  is  used  in 
structural  dynamics.  While  feasible  in  many  cases,  a  sequential  imposition  is  not  al¬ 
ways  possible,  leading  to  the  development  of  consistent  constraint  explicit  methods 
[14]  that  append  constraint  equations  to  the  mass  matrix  using  Lagrange  multipli¬ 
ers.  Although  this  method  is  effective,  it  is  also  more  expensive  in  comparison  to 
standard  explicit  formulations.  Constraints  that  cannot  be  analytically  eliminated 
or  imposed  sequentially  are  therefore  imposed  using  the  penalty  method.  If  a  con¬ 
straint  approach  is  taken  to  enforce  the  continuity  between  patches,  the  bending 
strip  is  therefore  preferred  based  on  efficiency  considerations. 


4  The  blended  shell  element  formulation 


The  blended  shell  formulation  proposed  here  may  be  thought  of  as  selectively 
adding  rotational  degrees  of  freedom  instead  of  adding  constraints.  Although  adding 
degrees  of  freedom  to  satisfy  a  set  of  constraints  sounds  counter  intuitive,  choos¬ 
ing  generalized  coordinates  that  automatically  satisfy  the  required  kinematic  con¬ 
straints  has  a  long  history  in  the  analysis  of  mechanisms  [29].  This  approach  elim¬ 
inates  the  need  to  enforce  the  nonlinear  C1  continuity  constraints,  eliminates  the 
possibility  of  locally  over  constraining  the  structure  at  complicated  intersections, 
and  simplifies  imposing  the  boundary  conditions.  More  specifically,  rotational  de¬ 
grees  of  freedom  are  selectively  added  to  a  rotation-free  formulation  [6]  as  needed 
to  allow  patches  to  be  connected.  This  formulation  is  not  restricted  to  the  control 
points  on  the  boundary,  but  permits  the  rotational  degrees  of  freedom  to  be  added 
anywhere  within  the  patch  to  permit  arbitrary  connections.  It  also  simplifies  the 
imposition  of  the  boundary  conditions. 

The  blended  formulation  may  be  interpreted  as  one  that  switches  from  a  Kirchhoff- 
Love  theory  to  a  Reissner-Mindlin  theory  at  appropriate  locations  in  the  structure. 
With  thin  shell  theories,  the  angle  between  two  patches  is  rigidly  fixed,  but  this 
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is  not  true  for  shear  deformable  theories.  The  ability  to  switch  between  theories 
therefore  has  some  attractive  properties  for  large  deformation  problems  such  as 
metal  forming.  One  operation  that  is  commonly  performed  on  metal  stamping  is 
“flanging,”  where  the  sheet  is  locally  folded  back  on  itself  to  stiffen  the  structure, 
introducing  a  6'°  discontinuity  along  the  fold  line.  This  operation  is  commonly 
simulated  in  metal  stamping  analyses  using  standard  C°  shell  elements  without  any 
difficulty,  but  similar  simulations  would  not  be  possible  with  a  pure  C 1  formulation. 


4.1  The  shell  kinematics 


The  shell  kinematics  follow  the  degenerated  solid  approach  to  shells  [17]  where  the 
volume  is  reduced  from  three  dimensions  to  a  two-dimensional  representation.  All 
coordinates  are  in  the  current  configuration.  By  assumption,  the  thin  dimension  is 
associated  with  the  third  parametric  coordinate  s:i.  For  convenience,  —  1  <  s3  < 
+1.  A  reference  surface  a;RS(si,  s2)  midway  between  the  upper  and  lower  surfaces 
has  the  parameteric  coordinate  s3  =  0.  A  fiber  vector  of  unit  length  f(s\,  s2)  and 
a  thickness  function  h(si,  s2)  are  introduced  so  that  the  coordinates  of  any  point  in 
the  shell  are  defined  in  terms  of  the  parametric  coordinates  as 

h 

*(si,s2,s3)  =  £CRS(si,s2)  +  s3-(s1,s2)f(s1,S2)  (8) 


In  most  shell  formulations,  the  thickness  function  is  assumed  to  be  independent  of 
time  and  therefore  volume  is  not  conserved  for  large  membrane  strains.  Recent  ef¬ 
forts  to  incorporate  thickness  strains  by  adding  extra  degrees  of  freedom  have  been 
very  successful,  e.g.,  [9,10].  Earlier  efforts  [16],  where  the  normal  strain  required  to 
satisfy  the  zero  normal  stress  condition  is  integrated  through  the  thickness  to  update 
the  thickness,  also  work  well.  For  notational  simplicity,  the  thickness  is  assumed  to 
be  independent  of  time,  and  the  extension  of  the  current  work  to  incorporate  thick¬ 
ness  changes  may  be  introduced  using  any  of  the  methods  developed  for  standard 
shell  elements. 


The  Kirchhoff-Love  hypothesis  defines  /  as  the  unit  normal  n  to  the  reference 
surface.  The  normal  is  defined  from  the  derivatives  of  the  midsurface, 


dx  dx 
P  =  ~ 

OSi  US2 


and  n 


P_ 
I P 


(9) 


With  this  assumption,  the  motion  of  the  shell  volume  is  uniquely  defined  by  the 
motion  of  the  reference  surface,  and  it,  in  turn,  is  defined  entirely  by  the  motion  of 
the  control  points  xa- 


The  Reissner-Mindlin  hypothesis  adds  the  coordinates  defining  the  orientation  of  f 
as  new  solution  variables,  allowing  the  development  of  transverse  shear  strain.  The 
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orientation  of  the  fiber  vector  can  be  represented  by  the  points  on  the  surface  of 
a  sphere,  giving  two  extra  rotational  degrees  of  freedom,  f{6 1,  02).  Working  with 
two  finite  angles  in  practice,  however,  is  difficult  in  structures  with  complicated 
shell  intersections.  Most  shell  element  formulations  use  the  simpler  representation 
of  the  three  degrees  of  freedom  associated  with  the  angular  velocity  oj  about  the 
global  coordinate  directions, 

f  =  ujxf.  (10) 

As  a  practical  matter,  many  Reissner-Mindlin  shell  elements  set  the  fiber  vector  to 
be  the  current  configuration  of  the  reference  surface  normal,  e.g.,  [3],  simplifying 
the  formulation  and  increasing  its  robustness. 


4.2  The  discrete  kinematics 


The  motion  of  the  reference  surface  (denoted  with  the  superscript  RS)  is  commonly 
expressed  as 

xrs  =  ^2Naxa,  ur*  =  J2naua,  (11) 

A  A 

where  uRS  is  the  displacement  of  the  reference  surface,  but  there  are  numerous 
choices  for  the  discrete  expression  of  the  displacement  contribution  associated  with 
the  fiber  vector.  Most  Kirchhoff-Love  elements  directly  evaluate  the  normal  using 
Equation  9  at  (si,  s2),  be.,  f(s i,  s2)  =  n(si,  s2)  in  Equation  8,  giving  the  discrete 
expression  for  the  coordinates, 

*(si,  s2,  s3)  =  na(si,  s2)xa  +  s3^(si,  s2)n(si,  s2).  (12) 

A  Z 

The  velocity  field  for  the  Kirchhoff-Love  formulation  is  obtained  by  directly  dif¬ 
ferentiating  Equation  12, 


x  =  J2NAxa  +  sj^-h 

A  1 

h  —  T — -  (I  —  n®n)p 
\p\ 

dxRS  dxRS 

(13) 

(14) 

v  =  „ —  x  — — 

ds i  ds2 

dxRS  dxRS  dxRS  dxRS 

(15) 

1  dsi  ''  d s2  '  dsi  ''  ds2 

(16) 

where  the  dependencies  on  the  parametric  coordinates  are  omitted  for  brevity. 
Using  the  degenerated  solid  approach  to  formulate  a  shell,  the  fiber  contribution  is 


10 


interpolated  from  the  nodes  as  in  the  Hughes-Liu  element  [17], 


f(sl,s2)  =  NA(si,s2)fA,  (17) 

A 

and  the  resulting  discrete  representation  of  the  coordinates  is 

*(sij  s2,  -S3)  =  na(si,  s2)  (xA  +  s3 .  (18) 

Our  previous  isogeometric  shells  [5,6]  use  Equation  18  with  f  A  =  nA  for  robust¬ 
ness  and  simplicity. 

Differentiating  Equation  18  for  the  velocity  field  using  the  normal  vector  for  the 
fiber  direction  gives 

X  =  J2na  [xa  +  •  (19) 

For  the  Reissner-Mindlin  isogeometric  shells  [5],  the  time  derivative  of  the  normal 
is  evaluated  as 

hA  =  u)Ax  nA,  (20) 

the  same  expression  used  in  standard  C°  elements.  The  time  derivative  of  the  nor¬ 
mal  may  also  be  evaluated  by  taking  the  time  derivative  of  Equation  9  evaluated  at 
A  as  in  the  rotation-free  formulation  [6], 


nA  = 


Pa  = 


Pa  = 


\Pa\ 

dx 

dsi 

dx 

dsi 


(■ I  ~nA®  nA)  pA 


\A  X 


\A  X 


dx 

ds2 

dx 

ds2 


dx 

dsi 


A  X 


dx 

ds2 


(21) 

(22) 

(23) 


Note  that  “evaluated  at  A”  means  that  for  every  control  point  A  there  is  a  unique 
location  in  the  parametric  domain  where  nA  is  evaluated  (see  [5,6]  for  further  de¬ 
tails). 


4.3  Selective  enrichment  with  the  blended  formulation 


The  only  difference  between  the  kinematics  of  the  rotation-free  formulation  [6]  and 
the  Reissner-Mindlin  formulation  [5]  is  the  evaluation  of  hA,  as  given  by  Equations 
21  and  20,  respectively.  Since  these  equations  refer  only  to  control  point  variables, 
the  equations  may  be  invoked  on  a  control  point  by  control  point  basis.  In  effect, 
this  shell  formulation  permits  the  smooth  blending  of  a  thin  shell  theory  with  a 
shear  deformable  one. 
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Rotational  degrees  of  freedom  may  therefore  be  selectively  introduced  to  permit 
the  transmission  of  the  bending  moment  along  lines  of  C°  continuity  without  intro¬ 
ducing  constraints  by  locally  using  a  shear  deformable  shell  theory.  This  is  accom¬ 
plished  by  using  the  kinematics  defined  by  Equations  18  and  19,  and  choosing  n  to 
be  defined  by  Equation  21  for  the  control  points  in  the  thin  shell  regions,  A  e  T, 
and  by  Equation  20  for  control  points  in  the  shear  deformable  regions,  d6(S, 


f  h  N 

x  =  J2na  (  xA  +  s3-nA 


x  =  J2  naXa 


h  (  1  ' 

■  s3  ~  I  51  Na] — |  (I  -nA®  nA)  pA  +  NaUJa  x  nA 

z  Vast  \Pa\  Aes  / 


(24) 


(25) 


Elements  with  control  points  belonging  to  both  sets  exhibit  a  behavior  somewhere 
between  the  pure  thin  shell  and  shear  deformable  shell  theories.  Example  problems 
presented  later  explore  the  influence  of  using  a  shear  deformable  theory  in  relatively 
localized  regions  of  an  otherwise  thin  shell  analysis. 


4.4  Implementation 


The  implementation  is  a  combination  of  the  Reissner-Mindlin  [5]  and  rotation-free 
[6]  isogeometric  shell  implementations.  These  shells,  in  turn,  follow  the  general¬ 
ized  element  strategy  [4].  The  focus  here  is  on  how  much  of  the  various  required 
calculations  can  be  carried  out  in  common  in  order  to  both  maximize  the  efficiency 
of  the  implementation  and  define  where  independent  operations  must  be  performed. 

While  the  primary  motivation  for  introducing  rotations  is  joining  patches  in  a  multi¬ 
patch  analysis,  they  are  also  convenient  for  imposing  boundary  conditions.  Addi¬ 
tionally,  there  are  occasions  when  having  C°  lines  running  through  a  patch  sim¬ 
plifies  the  model  generation.  To  incorporate  these  extensions,  the  definition  of  the 
set  S  is  generalized  to  include  any  control  point  with  rotational  degrees  of  freedom 
even  if  they  are  not  on  the  patch  boundary,  and  the  set  T  contains  the  remainder. 
Every  element  may  therefore  contain  an  arbitrary  combination  of  control  points 
with,  and  without,  rotational  degrees  of  freedom. 

Gauss  quadrature  is  performed  over  the  reference  surface  spanned  by  si  and  s2, 
and  through  the  thickness  direction  spanned  by  S3.  Many  of  the  common  terms 
are  not  functions  of  S3.  Efficiency  is  therefore  enhanced  by  evaluating  them  once 
before  performing  the  integration  loop  through  the  thickness,  and  making  the  in¬ 
tegration  over  the  reference  surface  the  outer  loop.  Variables  associated  with  the 
outer  and  inner  integration  loops  are  given  the  superscripts  G  and  g,  respectively. 
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The  integration  point  within  an  element  is,  for  example,  specified  by  the  parametric 
coordinates  (sf ,  sf ,  sf),  and  the  volume  integration  is  evaluated  as 


where  the  integration  weights  for  G  and  g  are  WG  and  w9,  respectively.  For  com¬ 
pactness,  the  evaluation  of  a  function  at  the  specific  integration  point  (sf ,  sf ,  sf)  is 
denoted  by  the  dual  superscripts  “Gg”  and  the  volume  contribution  of  integration 
point  Gg  is  denoted  V°9. 


4.5  Evaluation  of  the  normals  and  their  time  derivatives 


Since  the  basis  functions  are  not  interpolatory,  there  is  no  natural  location  for  eval¬ 
uating  the  normal.  The  derivatives  of  x  and  x  with  respect  to  the  two  parametric 
coordinates  si  and  s2  in  Equations  15  and  16  are  replaced  by  the  vectors  t,  and  i, 


tAi  =  CBiXB 
B 

tAi  =  CBiXB 
B 


(27) 

(28) 


where  G'r]t  are  evaluated  by  a  lifting  operation  (for  details  and  explicit  formulae, 
see  [6]).  The  normals  are  computed  at  each  control  point  outside  of  the  element 
integration  loop  since  they  are  independent  of  the  integration  points, 


Pa  —  t A\  X  f,42j  A 


Pa 

I  PaY 


no  sum  on  A. 


(29) 


The  time  derivative  of  the  normals  are  dependent  on  the  degrees  of  freedom  at  the 
control  point, 


tia  =  uja  x  nA,  A  e  S 

(30) 

nA=  (I  nAx  nA)pA,  AeT 

\Pa\ 

(31) 

PA  —  f  .41  X  t  A2  +  f  ,41  X  t  42 

(32) 

Evaluating  the  time  derivative  of  the  normal  once  outside  the  element  integration 
loops  simplifies  the  evaluation  of  the  velocity  gradient,  which  is  evaluated  at  every 
integration  point. 
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4.6  Evaluation  of  the  Jacobian  and  the  velocity  gradient 


The  summation  convention  is  used  as  follows:  When  repeated,  it  is  to  be  understood 
that  subscripts  i,  j,  k,  and  i  are  summed  over  1,  2,  and  3,  unless  it  is  explicitly 
indicated  to  sum  over  1  and  2. 

The  Jacobian  J  and  velocity  gradient  L  are  functions  of  the  control  point  coor¬ 
dinates,  the  normal,  and  the  time  derivative  of  the  normal,  and  they  are  therefore 
evaluated  without  regard  to  the  degrees  of  freedom  at  the  control  points.  Note  that 
both  are  affine  functions  of  S3,  allowing  most  of  the  calculations  to  be  performed 
once  at  each  reference  surface  integration  point. 


tG(  \  —  _  jGl 


+  s-vh 


G2 

ik 


J?*=Y.a-frx*  k  =  1’2 

JS‘=T.dfAf^t  k  =  1,2 


dsi-  2 


(33) 

(34) 

(35) 

(36) 


jG\  ,  „  rG2 


y 


+  s3  Li 


-1 


E  4, 

k= 1,2 

E  41 

k=  1,2 


E 


ONa. 
d  si 


%  Ai 


+  J 


~"Y.NAhfnAi 

.  A  Z 


3  3 


V-  dNA  hA  . 


(37) 

(38) 

(39) 


4. 7  Stress  evaluation 


The  stress  is  updated  in  a  co-rotational  formulation  similar  to  the  one  used  by  Bel- 
tyschko  and  Tsay  [3].  It  is  completely  independent  of  the  blending,  and  briefly 
documented  here  only  for  completeness.  A  local  co-rotational  coordinate  system 
ef,  i  =  1,  2,  3  is  defined  on  the  reference  surface  and  is  used  for  all  the  integration 
points  through  the  thickness.  The  normal  direction  defines  e|,  and  the  other  two  di¬ 
rections  are  calculated  using  the  procedure  of  Hughes  and  Liu  [17,15].  The  rotation 
matrix  from  the  local  to  the  global  coordinate  system,  v 9  =  Rve,  is 


R 


pi  pi 
C1  e2  e3 


(40) 
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The  rate  of  deformation  D  in  the  global  system  is  calculated  from  the  velocity 
gradient,  and  the  local  rate  De  is  obtained  by 

D°9  =  ^  ( Lg 9  +  ( LG9)t )  ,  D-  =  RTDGgR.  (41) 

The  normal  component  calculated  in  this  manner  will,  in  general,  not  result  in 
the  normal  stress  condition  ats  =  0  being  satisfied.  A  simultaneous  solution  for 
0-33  (.Df 3)  =  0  and  the  stress  update  is  performed  using  the  appropriate  algorithm 
[28,14,26].  Finally,  the  updated  stress  is  rotated  into  the  global  coordinate  system 
at  the  end  of  the  time  step, 

a°9  =  R(t(Rt.  (42) 


4.8  Evaluation  of  the  residual 


To  preserve  the  commonality  between  the  two  different  types  of  control  points  for 
as  long  as  possible,  the  virtual  velocity  is  expressed  in  terms  of  the  normal  and 
differentiated  to  obtain  the  expression  for  the  virtual  gradient, 


d5x  jdNA  hA  (  dNA  ds3\  \ 

^  Y  +  T  (S3lte  +  Njidtc ) SnA! 

Index  notation  simplifies  gathering  the  common  terms  in  the  sequel, 


(43) 


dSxi  srldNA\-  ,  llA  9Na  ,  at  dsAs- 

dXj  Y  l  dxi  2  V  &x3  dx3 


(44) 


+  hf4G%i 


where  J  1  is  the  inverse  of  the  Jacobian  matrix. 


Numerically  integrating  separately  in  the  reference  plane  si  —  s2  and  through  the 
thickness  direction  s3,  the  virtual  work  associated  with  the  stress  is 


SW  =  EE  E  4 Vc%l8ffsiA.v^ 

A  Gg  k=  1,2  USk 


E E 43(jGs)4 y  «^yGs 


(45) 


Separating  the  terms  based  on  their  dependencies  on  G  and  g,  the  virtual  work  is 


15 


flv  =  E»*E  E 

A  G  k=  1,2  9 

+Ey«A,{E  E  ^’Z4*%,(J°%1va’ 

A  z  G  k= 1,2  *  9 

+E^E^(^9)si1i/°!’}  (46) 

G  9 


Collecting  the  common  terms  that  are  summed  over  g  defines  two  resultant  vectors, 

*S  =  E'7.?VG%1'/G9  and  R?k  =  ^2s^(Jc%}Va‘.  (47) 

9  9 

Substituting  them  into  Equation  46  gives  the  forces  that  are  work  conjugate  to  5 x 
and  5h, 


=E  E  R 

G  k= 1,2 


GdNl 

ik  dsk 


hA 


,dN9 


f%  =  f  £(  Y,  ^k-i^  +  R%n9 


G  \fc=l,2 


<9st 


(48) 

(49) 


Note  that  the  residual  evaluation,  consisting  of  the  summation  over  the  two  nested 
loops  for  the  integration,  is  independent  of  whether  a  control  point  has  rotational 
degrees  of  freedom  or  not  to  this  point.  The  only  operations  that  depend  on  the 
degrees  of  freedom  are  performed  outside  of  the  integration  loop  just  before  the 
element  residual  assembly. 


For  the  rotation-free  control  points,  the  discrete  virtual  power  equation  at  control 
point  A  gives 

n  =  n+E^-^f-  (50) 

Introducing  the  projection  matrix  PAl]  and  the  common  terms  Cj]k  (see  Equation 
27),  the  final  expression  for  the  control  point  translational  forces  is 


Fll  =  Fii  +  Y.^niPBjkC^ 

B 


where 


R*Aij  — 


klAiJlAj  ) 


(51) 


(52) 


For  the  control  points  with  rotational  degrees  of  freedom,  the  discrete  virtual  power 
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equations  at  control  point  A  are 


Fm  = 


tax 
*  Ai 


and 


*M  = 


AjduAi  -e^nAkFAr 


p" 


(53) 


5  Example  calculations 


Explicit  finite  element  formulations  traditionally  rely  on  lower-order  elements  with 
uniformly  reduced  integration  for  their  computational  efficiency.  Therefore,  explicit 
analysis  typically  uses  elements  with  linear  basis  functions,  one  point  integration, 
and  hourglass  control.  Higher  order  Lagrange  elements  are  not  practical  for  explicit 
analysis  because  of  the  increased  cost  of  evaluating  the  element,  their  instability  un¬ 
der  large  deformations,  and  their  reduced  time  step  size  due  the  large  errors  in  their 
highest  frequencies.  The  convergence  of  NURBS  basis  functions  at  their  highest 
frequencies  eliminates  the  time  step  size  penalty,  but  the  cost  of  evaluating  the  ele¬ 
ments  still  increases  with  p.  Based  on  previous  benchmarks  [4],  quadratic  NURBS 
elements  with  uniform  reduced  integration  offer  the  greatest  opportunity  for  en¬ 
hancing  the  spatial  accuracy  without  increasing  the  analysis  cost  in  comparison  to 
current  formulations.  Throughout,  each  element  is  integrated  with  the  2x2  uni¬ 
formly  reduced  Gauss  quadrature  rule  in  the  plane,  and  the  three-point  rule  through 
the  thickness.  Indeed,  as  shown  previously  [4]  and  in  the  last  example  in  this  sec¬ 
tion,  isogeometric  analysis  with  quadratic  NURBS  offers  the  possibility  of  increas¬ 
ing  accuracy  at  a  reduced  computational  cost.  The  focus  in  the  examples  presented 
here  is  therefore  on  quadratic  NURBS  exclusively.  However,  we  note  that  recent 
work  on  isogeometric  collocation  methods  in  explicit  dynamics  suggests  that  this 
picture  may  change. 


5.1  Nonlinear  L-b earn 


In  order  to  demonstrate  the  ability  of  the  current  formulation  to  preserve  the  initial 
angle  of  multiple  patches  along  C°  lines,  an  L-shaped  cantilever  beam  subjected  to 
a  point  load  is  chosen.  This  example  has  been  analyzed  by  Kiendl  et  al.  [20,21]  to 
show  the  performance  of  the  constraint  equations  and  the  bending  strip  method  in 
multi-patch  analysis. 

Quadratic  NURBS  elements  in  three  different  model  configurations  are  analyzed 
as  shown  in  Figure  2.  Starting  from  the  left,  the  first  model  has  two  patches  with 
rotation-free  isogeometric  shell  elements  everywhere.  The  second  model  is  split 
into  two  patches  with  the  blended  isogeometric  shell  formulation  with  the  rota¬ 
tional  degrees  of  freedom  added  along  the  C°  edge.  Finally,  a  two  patch  model 
with  the  Reissner-Mindlin  isogeometric  shell  formulation  [5]  is  also  analyzed.  The 
structural  response  is  shown  in  Figure  3.  The  first  model,  as  expected,  illustrates  the 
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Fig.  2.  Initial  meshes  of  the  L-beams  modeled  with  quadratic  NURBS.  From  left  to  right: 
two  patches  of  rotation-free  elements,  two  patches  of  blended  elements  with  rotations  along 
the  fold  line,  and  finally  two  patches  of  Reissner-Mindlin  shell  elements. 
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Fig.  3.  Final  configuration  of  the  L-beams.  From  left  to  right:  two  patches  of  rotation-free 
elements,  two  patches  of  blended  elements  with  rotations  along  the  fold  line,  and  finally 
two  patches  of  Reissner-Mindlin  shell  elements. 

inability  of  rotation-free  shells  to  transmit  bending  moments  across  C°  lines.  The 
C°  line  behaves  like  a  hinge  and  the  initial  90°  angle  is  not  maintained.  The  other 
models  give  identical  results,  demonstrating  the  proper  transmission  of  bending 
moments  with  the  blended  formulation. 


5.2  Vibration  of  a  thin  simply-supported  square  plate 


The  previous  example  evaluated  the  accuracy  of  different  methods  for  multi-patch 
analysis  subject  to  single  modes  of  deformation.  An  eigenvalue  analysis  introduces 
a  large  number  of  modes  of  deformation  with  the  complexity  of  the  deformations 
increasing  with  the  frequency,  giving  additional  insight.  This  example  was  used  to 
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Fig.  4.  The  geometry  of  the  local  constraint  equation.  The  C°  interface  is  normal  to  the 
plane  of  the  figure  and  located  at  the  middle  control  point.  Aw  is  in  the  global  coordinate 
system  perpendicular  to  the  plate,  and  Xe  is  in  the  plane  of  the  plate  and  perpendicular  to 
the  C°  boundary. 

test  isogeometric  analysis  in  [13],  and  was  also  analyzed  in  [6]  to  demonstrate  the 
accuracy  of  the  rotation-free  isogeometric  shell  formulation  in  linear  analysis.  The 
exact  eigenvalues  in  radians  for  a  square  plate  of  length  L  and  thickness  h,  using 
thin  plate  theory,  are 

Wij  =  C(i2+j2)  0  <i,j 


p  (12  (1  —  z/2))  L2 

where  E  is  Young’s  modulus,  v  is  Poissons’s  ratio,  and  p  is  the  density.  In  this 
example,  the  values  are  chosen  to  be  E  =  107,  v  =  0.3,  p  =  1.0,  h  =  0.05, 
L  =  10. 


(54) 

(55) 


The  square  plate  is  evenly  divided  into  two  patches  and  the  eigenvalue  problem  is 
solved  using  three  formulations: 

(1)  The  rotation-free-formulation  [6]  with  the  geometric  continuity  enforced  with 
linear  constraints  formulated  in  a  manner  similar  to  Kiendl  et  al.  [21].  The 
constraints  are  discussed  in  more  detail  below.  This  formulation  is  referred  to 
as  RF1. 

(2)  A  second  rotation-free  formulation,  based  on  the  updated  Lagrangian  formu¬ 
lation  given  by  Equation  12  with  C1  continuity  enforced  using  the  same  linear 
constraints  as  in  (1).  The  two  rotation-free  formulations  differ  in  where  their 
normals  are  evaluated,  however  they  give  the  same  eigenvalues  on  a  single 
patch  for  the  square  plate  problem.  Of  interest  here  is  the  sensitivity  of  the 
eigenvalue  accuracy  to  the  different  rotation-free  formulations  when  using  the 
constraints  to  enforce  continuity.  This  formulation  is  referred  to  as  RF2. 

(3)  The  blended  formulation  with  the  common  boundaries  between  the  patches 
having  translational  and  rotational  degrees  of  freedom  while  the  interior  of 
the  patches  is  rotation-free.  Since  constraints  on  the  rotations  are  not  imposed 
on  the  perimeter  of  the  plate,  the  perimeter  control  points  have  only  the  three 
constrained  translational  degrees  of  freedom. 

The  linear  constraints  enforcing  the  continuity  [21]  are  defined  in  terms  of  the  local 
reference  geometry  and  the  infinitesimal  displacements  as  illustrated  in  Figure  4. 
Neglecting  the  higher  order  terms,  the  linear  constraint  equation  in  the  local  co- 
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Fig.  5.  Convergence  of  the  first  eigenvalue  for  the  two  patch  mesh  of  the  square  plate 
problem  as  a  function  of  the  number  of  control  points. 


Fig.  6.  Convergence  of  the  first  eigenvalue  for  the  four  patch  mesh  of  the  square  plate 
problem  as  a  function  of  the  number  of  control  points. 

ordinate  system  defined  at  the  shared  boundary  of  two  coplanar  patches  simplifies 
to 

AXfAw2  -  AX^AWl  =  0.  (56) 

where  Xe  is  in  the  local  coordinate  system  of  the  flat  plate  with  one  axis  aligned 
with  the  patch  boundary. 

The  convergence  rates  for  the  first  eigenvalue,  using  two  quadratic  NURBS  patches, 
are  shown  in  Figure  5. 

An  additional  C°  discontinuity  is  introduced  by  splitting  the  domain  into  four 
square  patches.  This  problem  is  fundamentally  different  from  the  previous  one  be¬ 
cause  it  has  intersecting  C°  fines  involving  nine  control  points  in  the  constraints  at 
the  center  of  the  plate.  If  the  constraint  patterns  along  the  two  C°  fines  are  blindly 
imposed  along  the  C°  fines,  the  nine  control  points  are  subjected  to  six  constraint 
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equations,  producing  a  singular  system.  Deleting  both  constraints  for  the  center 
control  point  improved  the  accuracy  of  some  modes  and  made  the  accuracy  of  oth¬ 
ers  worse  compared  to  the  single  patch  results,  but  deleting  only  one  restored  the 
accuracy  to  almost  exactly  the  single  patch  results.  Constraints  making  the  center 
control  point  displacement  the  average  of  a  set  of  the  surrounding  control  points, 

WC  -  T7  WA  =  °,  (57) 

Aec 

where  C  designates  the  central  control  point  and  C  is  a  subset  of  the  control  points 
surrounding  C,  worked  equally  well.  General,  accurate,  robust  strategies  for  con¬ 
straining  control  points  at  arbitrary  shell  intersections  are  currently  unavailable. 

No  special  treatment  was  required  for  the  blended  formulation.  The  center  node  was 
treated  like  any  other  node  at  a  C°  line,  i.e.,  it  had  a  full  set  of  rotational  degrees  of 
freedom.  The  control  points  on  the  simply  supported  boundary  again  had  only  the 
three  constrained  translational  degrees  of  freedom. 

The  convergence  rates  are  shown  in  Figure  6.  These  results  are  virtually  identical 
to  the  results  for  one  and  two  patches. 


5.3  Gravity-loaded  L- shaped  plate 


An  L- shaped  plate  is  loaded  by  gravity.  The  problem  is  specified  numerically  and 
the  geometry,  material,  and  boundary  condition  data  are  provided  in  Figure  8.  The 
parameters  are  chosen  to  make  this  a  nonlinear,  finite  deformation  problem  with 
coupled  membrane  and  bending  behavior.  The  plate  is  discretized  with  NURBS, 
and  two  meshes,  coarse  and  fine,  are  shown  in  Figure  9.  The  parameterization  has 
a  discontinuous  derivative  along  the  line  of  symmetry.  As  a  result,  the  underlying 
discretization  is  only  G'°-continuous  along  the  line  of  symmetry  and  a  direct  appli¬ 
cation  of  a  rotation-free  shell  formulation  is  not  possible  in  this  case. 

One  approach  to  overcome  this  difficulty  is  to  use  the  bending  strip  method.  A 
bending  strip  patch  may  be  generated  along  the  symmetry  line.  A  typical  element 
of  that  patch  is  depicted  in  Figure  10.  The  figure  reveals  the  difficulty  with  the  ap¬ 
plication  of  the  bending  strip  method  in  this  case:  Due  to  the  abrupt  change  in  the 
parameterization  near  the  symmetry  line,  the  covariant  basis  vector  Gi,  which  is 
typically  used  to  define  the  direction  of  the  bending  stiffness  operator,  is  rapidly 
changing  direction  from  one  side  of  the  bending  strip  to  the  other.  Furthermore, 
almost  everyhere  Gi  has  a  component  parallel  to  the  interface,  which  will  incor¬ 
rectly  stiffen  the  structure  in  this  direction.  Despite  this  difficulty,  the  bending  strip 
method  was  made  to  work  in  this  case  by  integrating  the  contribution  of  the  bending 
strip  terms  with  one-point  quadrature  placed  in  the  middle  of  the  element.  At  this 
quadrature  point  the  covariant  basis  vectors  are  orthogonal,  and  Gi  points  in  the  di- 
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All  constraints  (locks) 


Two  constraints  removed 


One  constraint  removed 


Mean  constraint  pattern  T1  Mean  constraint  pattern  T2Mean  constraint  pattern  T3 


Fig.  7.  Constraint  configurations  for  the  four-patch  mesh.  Solid  lines  indicate  constraints 
of  the  type  defined  by  Equation  56  and  dashed  lines,  Equation  57.  Control  points  shared 
between  adjacent  patches  are  indicated  with  white  dots,  and  those  associated  with  only  one 
patch,  black  dots.  The  configuration  with  all  constraints  locks.  With  only  one  constraint 
removed,  the  model  works  correctly.  The  three  patterns  is  the  second  row  using  Equation 
57  worked  equally  well. 


Fig.  8.  L-shaped  plate.  Problem  setup  and  data. 

rection  orthogonal  to  the  interface.  Figure  1 1  shows  the  results  of  the  computations 
with  and  without  the  bending  strip. 


This  “fix”  is  somewhat  ad-hoc,  and,  in  what  follows,  the  L-shaped  plate  example 
is  used  to  illustrate  the  good  behavior  of  the  proposed  blended  formulation,  which 


Fig.  9.  L-shaped  plate.  NURBS  meshes.  Note  that  the  parameterization  of  the  domain  has 
a  discontinuous  derivative  along  the  symmetry  line. 


Fig.  10.  L-shaped  plate.  Illustration  of  the  difficulty  with  the  application  of  the  bending  strip 
method  to  this  example.  The  direction  of  the  covariant  basis  vector  Gi,  which  is  typically 
taken  as  the  direction  for  the  bending  stiffness  operator,  is  only  truly  orthogonal  to  the 
interface  at  the  midpoint  of  the  parametric  line.  Also  note  that  the  covariant  basis  vectors 
are  orthogonal  only  at  this  location.  The  six  control  point  element  in  a)  overlays  the  two 
yellow  elements  in  b). 

does  not  rely  on  a  good  parameterization  near  C'°-continuous  interfaces. 

An  analytical  solution  is  not  available  for  this  problem.  In  a  one-patch  convergence 
study  with  p  —  4,  the  maximum  displacement  converged  to  0.003277,  the  value 
we  used  in  lieu  of  an  exact  solution  in  Figure  12.  The  rotation-free  solutions  (the 
bending  strip  method  uses  a  rotation-free  shell  [20])  converge  from  the  soft  side, 
while  those  with  rotational  degrees  of  freedom  converge  from  the  stiff  side. 

To  evaluate  the  benefits  of  the  blended  formulation  for  implicit  problems,  the  prob¬ 
lem  was  also  solved  using  the  Reissner-Mindlin  formulation  [5]  on  the  finest  mesh, 
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(a)  (b) 


Fig.  11.  L-shaped  plate.  Deformed  configuration  (scale  factor  of  200  is  employed)  colored 
by  the  z— displacement.  (a)  Result  without  the  bending  strip,  (b)  Result  with  the  bend¬ 
ing  strip  integrated  using  one-point  quadrature.  The  computation  without  the  bending  strip 
gives  an  unphysical  “kink”  at  the  interface,  while  the  one-point-quadrature  bending  strip 
formulation  gives  a  physically  correct  answer. 


Fig.  12.  Convergence  of  the  z  displacement  at  the  origin  for  the  L-shaped  domain. 

and  the  results  are  summarized  in  Table  1 .  The  reported  CPU  times  were  obtained 
with  the  standard  Linux  timing  function,  and  are  accurate  to  approximately  0.01 
second.  To  solve  the  nonlinear  problem  required  one  stiffness  factorization  and  six 
solves.  The  storage  is  as  expected:  The  blended  formulation  has  half  the  number 
of  degrees  of  freedom  and  requires  one  quarter  the  storage.  For  a  flat  plate,  the 
Reissner-Mindlin  formulation  is  effectively  a  five  degree  of  freedom  formulation, 
with  the  normal  rotations  assigned  a  small  values  in  the  stiffness  matrix  to  avoid  be¬ 
ing  singular.  Although  the  concept  of  a  “bandwidth”  is  not  appropriate  for  a  sparse 
solver,  the  logically  regular  structure  of  the  mesh  leads  to  an  ordering  such  that 
the  ratio  of  the  Reissner-Mindlin  to  the  blended  formulation  factorization  costs  is 
very  close  to  (5/3)3  =  4.63,  and  for  the  solves,  the  ratio  is  roughly  (5/3)2  =  2.78. 
Although  the  Reissner-Mindlin  formulation  has  simpler  terms  for  the  stiffness  ma- 
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Table  1 

Comparison  of  Implicit  Reissner-Mindlin  and  Blended  Formulation  Solution  Costs 


Reissner-Mindlin 

Blended 

Ratio  (R-M)/B 

Number  of  equations 

51093 

25353 

2.01 

Memory  (103  8-byte  words) 

3735 

937 

3.98 

Element  CPU  (sec.) 

7.37 

5.65 

1.30 

Assembly  CPU  (sec.) 

1.29 

1.09 

1.18 

Factorization  CPU  (sec.) 

12.7 

2.74 

4.63 

Total  Solve  CPU  (sec.) 

0.27 

0.09 

3.00 

Total  CPU(sec.) 

21.63 

9.57 

2.26 

L 

l 

hp 

h 

E 

v 

P 

ay 

Eh 

V 


320  mm 
70  mm 
67.5  mm 
1.2  mm 
1.994  GPa 
0.3 

7850  kg/m3 
3.366  x  102  MPa 
1.0  MPa 
5.646  m/s 


Fig.  13.  Buckling  of  a  square  tube:  problem  description  and  the  mesh.  Only  one  quarter  of 
the  geometry  is  modeled  with  appropriate  symmetry  boundary  conditions.  Reprinted  from 
[5] 


trix  and  the  residual,  having  more  of  them  puts  it  at  a  disadvantage  in  terms  of  the 
overall  element  cost.  The  total  gain  in  speed  is  over  a  factor  of  two  for  this  modest 
sized  problem. 


25 


Fig.  14.  Final  configuration  of  the  square  tube  buckling  problem  using  a)  the  Belytschko-T- 
say  element,  b)  the  fully-integrated  quadrilateral  shell  element,  c)  the  p  =  2  isogeometric 
Reissner-Mindlin  shell,  d)  the  p  =  2  rotation-free  isogeometric  shell  with  the  continuity 
across  the  C°  lines  enforced  with  rigid  bodies,  and  e)  the  blended  element. 

Table  2 

Square  tube  buckling  benchmark  summary 


Element 

Type 

Number  of 

Elements 

Number  of 

Time  Steps 

CPU 

Seconds 

1 -Point  Quad 

2560 

352631 

833 

Full  Int.  Quad 

2560 

414148 

2861 

Isogeo.  RM 

640 

153596 

582 

Isogeo.  RF 

640 

173706 

835 

Isogeo.  blended 

640 

158192 

754 

5.4  Buckling  of  a  square  tube 


The  problem  of  accordion-mode  buckling  of  a  square  tube  [8]  was  analyzed  previ¬ 
ously  [5]  to  evaluate  quadratic  and  quartic  NURBS-based  isogeometric  Reissner- 
Mindlin  shell  elements.  The  problem  definition,  geometry,  material  parameters, 
and  the  NURBS  mesh  for  one  quarter  of  the  domain  are  shown  in  Figure  13.  An 
isotropic  elastic-plastic  material  with  linear  plastic  hardening  is  used  to  model  the 
material  response.  The  deformation  is  driven  with  a  constant  velocity  at  one  end  of 
the  tube  with  the  other  end  fixed.  A  geometric  imperfection  with  an  amplitude  of 
0.05  mm  triggers  the  buckling  at  a  height  of  67.5  mm  from  the  base. 

One  quarter  of  the  tube  is  modeled  using  two  quadratic  NURBS  patches  that  share 
the  control  points  along  the  common  edge,  leading  to  640  elements.  Appropriate 
symmetry  boundary  conditions  are  used  together  with  a  standard  single  surface 
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contact  algorithm  in  LS-DYNA  [14,8]  applied  to  the  automatically  generated  in¬ 
terpolation  elements  [4].  Each  quadratic  NURBS  element  is  subdivided  into  2x2 
interpolation  elements. 


Three  different  NURBS-based  isogeometric  shell  formulations  were  analyzed:  The 
Reissner-Mindlin  shell  [5],  the  rotation-free  shell  [6],  and  the  blended  shell.  Conti¬ 
nuity  for  the  rotation-free  formulation  is  enforced  between  the  patches  using  rigid 
bodies  [7]  composed  of  the  common  control  point  between  adjacent  patches  and 
the  first  interior  control  points  for  each  patch,  and  in  a  similar  manner,  their  rota¬ 
tional  boundary  conditions  at  the  upper  and  lower  boundaries  were  imposed  with 
rigid  bodies  composed  of  the  boundary  nodes  and  the  first  rows  of  nodes  interior  to 
the  mesh.  For  the  blended  formulation,  the  rotational  degrees  of  freedom  are  placed 
along  all  boundaries  of  the  two  NURBS  patches. 


The  final  deformations  for  the  three  shell  formulations,  and  the  LS-DYNA  solutions 
using  the  Belytschko-Tsay  element  [3]  and  the  type  16  fully-integrated,  4-node 
shell,  are  shown  in  Figure  14.  The  meshes  for  the  standard  elements  have  2560  ele¬ 
ments,  four  times  as  many  as  the  quadratic  NURBS.  Using  reduced  2x2  integration 
with  the  quadratic  NURBS  elements  results  in  the  three  isogeometric  models  hav¬ 
ing  the  same  number  of  integration  points  as  the  Belytschko-Tsay  model,  while  the 
fully-integrated  shell  has  four  times  as  many.  The  final  configurations  are  remark¬ 
ably  similar  given  the  five  different  formulations. 


Since  the  membrane  strains  are  small,  the  use  of  the  rigid  bodies  to  impose  the  con¬ 
straints  on  the  rotation-free  elements  does  not  introduce  a  noticeable  error.  Details 
of  the  rigid  body  formulation  are  given  in  the  Appendix.  The  rigid  bodies  accounted 
for  approximately  three  percent  to  the  total  CPU  time,  indicating  that  they  are  an  ef¬ 
ficient  approach  for  enforcing  continuity  between  patches.  The  blended  solution  is 
less  costly  than  the  rotation-free  solution  primarily  because  of  the  reduced  number 
of  time  steps.  On  an  element  per  time  step  basis,  the  cost  of  the  two  formulations  is 
nearly  identical. 


Table  2  summarizes  the  problem  sizes  and  the  CPU  times  for  a  single  Xeon  proces¬ 
sor  running  the  analyses  in  double  precision.  Comparing  the  computational  costs  of 
the  isogeometric  elements  to  the  traditional  elements,  it  is  clear  that  the  quadratic 
isogeometric  elements  are  less  costly  primarily  because  of  their  larger  stable  time 
step  size.  The  fully-integrated  quadrilateral  is  the  most  expensive  solution  because 
it  has  a  smaller  time  step  size  and  four  times  as  many  integration  points  in  the  mesh 
as  the  quadratic  NURBS.  In  this  calculation,  the  Reissner-Mindlin  formulation  is 
the  most  accurate.  However,  we  note  that  in  implicit  calculations  there  should  be 
significant  efficiency  advantages  in  the  blended  formulation. 
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6  Summary  and  conclusions 


There  are  significant  potential  savings  in  computational  cost  with  rotation-free 
shells  in  implicit  analysis.  Unfortunately,  actually  achieving  such  savings  requires 
dealing  with  required  constraints  at  shell  intersections,  along  folds  and  boundaries, 
and  at  the  interfaces  of  NURBS  patches.  No  general  strategy  exists  for  implement¬ 
ing  the  variety  of  situations  that  occur  in  practical  engineering  calculations.  On 
the  other  hand,  the  classical  Reissner-Mindlin  element  formulation  does  not  incur 
any  problems  of  this  kind,  at  the  price  of  twice  as  many  degrees  of  freedom.  In 
this  work  we  have  proposed  a  new  shell  formulation  that  blends  thin-plate  rotation- 
free  theory  with  Reissner-Mindlin  theory.  By  only  selecting  the  Reissner-Mindlin 
formulation  locally,  where  the  constraints  need  to  be  enforced,  we  produce  a  for¬ 
mulation  that  yields  the  savings  of  the  rotation-free  formulation,  combined  with  the 
robustness  and  generality  of  the  Reissner-Mindlin  formulation.  We  have  tested  the 
formulation  on  several  problems  and  it  has  performed  well  in  all  cases.  We  have 
also  identified  some  of  the  difficulties  encountered  implementing  constraint  equa¬ 
tions  in  rotation-free  calculations.  Our  computations  have  revealed  that  uniformly 
reduced  integration  quadratic  NURBS  elements  are  computationally  efficient  and, 
for  the  same  level  of  accuracy,  they  rival  the  fastest  known  low-order  one-point 
quadrature  shell  elements  in  speed.  We  believe  that  the  isogeometric  blended  shell 
formulation  has  the  potential  to  become  the  method  of  choice  for  a  wide  variety  of 
practical  engineering  shell  problems. 
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A  A  Simple  Rigid  Body  Dynamics  Formulation  for  Explicit  Calculations 


Rigid  bodies  are  a  standard  feature  in  LS-DYNA  [14].  The  algorithm  [14,7]  for 
updating  their  motion  is  summarized  here. 

The  equations  of  motion  for  the  center  of  mass  for  a  rigid  body  are 

m'xcM  =  F  and  Jib  =  M  —  oj  x  Jiu  (A.l) 

where  m  is  the  mass  of  the  body,  J  is  the  mass  moment  of  inertia  tensor,  xcm  are 
the  coordinates  of  the  center  of  mass,  lo  is  the  angular  velocity,  and  F  and  M  are 
the  sums  of  the  externally  applied  forces  and  moments,  respectively. 

The  mass  and  the  initial  center  of  mass  coordinates  are 


m  =  J2  mA  (A. 2) 

A 

xCm  =  —  Y\  Maxa  (A. 3) 

mA 

where  the  summation  is  performed  over  all  the  control  points  (or  nodes)  in  the  rigid 
body.  After  determining  the  initial  center  of  mass,  the  mass  moment  of  inertia  J  is 
calculated  as 

J  =  J2'mA(\\x a  -  xCM\\ 2 1  -  ( xA  -  xCM )  ®  ( xA  ~  xCM ))  +  JAI-  (A.4) 

A 

At  tn,  the  nodal  forces  F'\  and  moments  M\  are  summed  to  calculate  the  rigid 
body  forces  and  moments, 


Fn  =  J2Fa  (A. 5) 

A 

M"  =  EW- *cm)  X  n  +  (A.6) 


A 

and  they  are  substituted  into  Equation  A.l. 

Central  difference  time  integration  updates  the  rigid  body  velocity  and  the  center 
of  mass  coordinates. 


.  n+l/2  _  ■  n- 1/2  a  j.n  -  n 

xcm  —  xcm  a-  lal  xcm 

a,«+V2  =  ^n-1/2  +  Atn^n 


rj.n+1  _  vn 

xcm — xcm 


i  Aj.n+l/2  •  n+l/2 
A-  AAl  XCM 


(A. 7) 
(A. 8) 
(A. 9) 
(A.  10) 
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The  calculation  of  the  incremental  rotation  vector 


A6»n+1  =  A  tn+1/2un+1/2  (A.ll) 

is  preferred  to  integrating  a  finite  rotation  measure  for  simplicity  and  robustness. 
Traditional  finite  rotation  measures  consisting  of  three  angles  have  singular  posi¬ 
tions  and  quaternions  require  the  imposition  of  a  quadratic  constraint  [29]. The  in¬ 
cremental  rotation  matrix  A R  is  calculated  from  A 9n+1  with  the  Hughes-Winget 
formula  [18]  and  used  to  incrementally  rotate  the  mass  moment  of  inertia  tensor 
from  n  to  n  +  1. 

Finally,  the  control  point  velocity  and  coordinates  are  updated  to  n  + 1/2  and  n  + 1, 
respectively,  in  a  manner  guaranteed  to  exactly  preserve  the  rigid  body  behavior  of 
the  system. 


•  n+l/2 


X 


X 


n+1 

A 


(A. 12) 

xnA  +  Atn+1/2xA+1/2 . 

(A. 13) 
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